NASA Technical Memorandum 89059 


ESTIMATION OF BIAS ERRORS IN MEASURED 
AIRPLANE RESPONSES USING MAXIMUM 
LIKELIHOOD METHOD 

(NASA-TM-89059) ESTIMATION 01 HAS ERRORS 
IN MEASURED AIRPLANE RESPONSES USING MAXIMUM 
I IKELIHOCD METHOD (NASA) 34 p CSCL OlC 


G3/08 


Vladislav Klein 
Dan R. Morgan 


January 1987 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 


N87-15S77 

Unclas 
43 367 


SUMMARY 


A procedure for compatibility check of measured airplane responses is 
presented. This procedure includes estimation of bias errors in the measured 
data in terms of constant measurement biases and scale factors, and a compari- 
son of reconstructed responses with those measured. The model relating air- 
plane states and outputs is based on six-degree-of -freedom kinematic 
equations. In these equations the input variables are replaced by their 
measured values which are assumed to be without random errors. 

A maximum likelihood method is used as the estimation technique. The 
resulting algorithm is verified with simulated data and data from flight 
testing. The results from simulated data show that the increased number of 
unknown parameters and the correlation among them can degrade the accuracy of 
the estimates; however, moderate measurement noise level in the input vari- 
ables has only a small effect on the estimates. The maximum likelihood esti- 
mates from flight data were compared with those obtained by using an extended 
Kalman filter and a nonlinear fixed-interval smoother. This comparison showed 
no major differences in results of all three techniques. 
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Sy, a z longitudinal, lateral, and vertical accelerations, m/sec 
constant bias error in variable y 
expected value 
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acceleration due to gravity, m/sec 
altitude, m 

nonlinear output vector used to represent measurement 
sys tem 

cost function 
information matrix 

modified information matrix, see Eq. (17) 
element of inverse information matrix 
number of data points 
measurement - noise vector 
number of measured output variables 
number of unknown parameters 
p, q, r roll, pitch, and yaw velocities, rad/sec or deg/sec 

R measurement - noise covariance matrix 

s( ) standard error estimate 

t time 

u, v, w longitudinal, lateral, and vertical airspeed components, 
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V true airspeed, m/sec 
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x b' ^b' z b linear position coordinates of aircraft, m 
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position coordinates of a wind vane with respect to 
aircraft center of gravity, m 

position coordinates of 0 wind vane with respect to 
aircraft center of gravity, m 

output vector 

transformation matrix 

measurement vector 

angle of attack, rad or deg 

sideslip angle, rad or deg 

Knonecker delta 

convergence criteria, see Eq. (16) 
input vector 

vector of unknown parameters 

pitch angle, rad or deg 

scale factor error of variable y 

residual vector 

process - noise vector 

variance 

roll angle, rad or deg 


Subscripts : 


b 

E 


N 


0 

R 

v 


body axes 
measured quantity 
nominal value 
initial value 
uncorrected for bias error 
wind vane 


Matrix Exponents : 


T transpose matrix 

-1 inverse matrix 

Mathematical Notation: 

• over symbols denotes derivative with respect to time 

* over symbol denotes estimated value 

A incremental value 

Abbreviations : 


EKF 

ML 

NFIS 


extended Kalman filter 

maximum likelihood 

nonlinear fixed-interval smoother 
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INTRODUCTION 


For more than fifteen years there have been numerous attempts to estimate 
airplane states from measured flight data. This estimation is possible 
because of well known kinematic equations relating the airplane states and 
output variables. In many cases the estimated states have been used for data 
compatibility checks, that is, for comparison of measured and predicted 
response variables of an airplane. Because the measured data are corrupted by 
random and systematic errors, it was recognized that the state estimates 
should be combined with the estimation of unknown biases (parameters) to 
obtain satisfactory results. Various methods for state and parameter esti- 
mation were applied. They can be divided into two groups: 

1 . Methods for estimating separately unknown parameters and states 
(refs. 1 to 6) . For parameter estimation the maximum likelihood or nonlinear 
least squares techniques are used. For state estimation the system is assumed 
to be deterministic. That is, the state variables are simply obtained by the 
integration of model equations. 

2. Methods for estimating states and parameters simultaneously using an 
extended Kalman filter (refs. 7 and 10), or nonlinear smoother (ref. 11). A 
review of various approaches to the problem of airplane state estimation is 
presented in reference 12. 

The purpose of this report is (1) to develop a maximum likelihood algor- 
ithm applicable to general motion of an airplane; (2) to compile an efficient 
computer program based on this algorithm; and (3) to verify both the algorithm 
and program on simulated and real flight data. The report starts with the 
formulation of model equations and estimation techniques. Then several 
examples are presented. When the real flight data are analyzed the maximum 
likelihood results are compared with those obtained by using an extended 
Kalman filter and nonlinear smoother. 

MODEL EQUATIONS 

The mathematical model used for the data compatibility check is described 
by three sets of kinematic equations with the state variables consisting of 
three linear velocities u, v, and w; three Euler angles <(>, 9, and \(i ; and 
three linear positions x^, y^, and z^. The input variables in these equations 

are the linear accelerations a x , a , and a z and angular velocities p, q, 

and r. The form of the kinematic (state) equations considered can be found in 
various references (see, e.g. reference 9) . 

The following variables are measured: 

1. The inputs to the system a v , a , a , p, q, and r. 

x y t. 

2. The airspeed V, two incidence angles 9 and a 
angles <J>,0 , and \p and altitude h= — z 


, three Euler 



These variables represent the output of the system. The measured vari- 
ables z are corrupted by systematic and random errors. It is assumed that 
each of them can be expressed as 

z=(1+X)y+b+n (1) 

y y y 


where y is the true value of the output, 
by is the constant bias error, and n y is 


A is the unknown scale factor error, 

y 

the measurement noise. It is further 


assumed that the scale factor error is equal to zero for all the input vari- 
ables. This assumption will simplify an estimation procedure for remaining 
scale factor and hias errors. 


The system of state equations is simplified by deleting the equations for 
x^ and y^. Then, replacing the input variables in the remaining state equa- 
tions by their measured values results in the following set of state equations 
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The output equations take the form 


V_ = ( 1 + A ) /u 2 + v 2 + w 2 + b 
R v v 


8* = (1 + A Jtan 1 [ 

vR 8 L u - r 


V + r R, E x 8~ P R,E Z 8 
R,E y 8 + q R,E Z 8 




_ -1 r V + r R,E x 8 + p R,E z 8i 

8 = (1 + A.) tan cos a _ _+ b 

vR 8 L u - r R ^y 3+ vR,E 
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-1 r W " VeV P R,E y a 
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h R= (1 + V h 
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V (1 + V 9 + b e 
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(3) 


In Eq. (3) 8 and a v are the incidence angles measured by the vane. The 
sideslip angle 0* is measured in the Ox^y^ plane, whereas 0 is measured 
between the wind vector and its projecton on the Ox^^ plane. These two 
angles are therefore defined as 

-1 v 
0* = tan — 
u 

„ . -1 v 

8 = sin - 


Sometimes the measured sideslip angle has been corrected for the c.g. offset 
prior to the compatibility check. In such case 3 can be computed from the 

simplified equation as indicated above. In the output equations for h^ 

and i[i R the constant bias terms are omitted because the reference values 
for \J> and h can be selected arbitrarily. Further, in Egs . (2) and (3), R 
indicates the variable uncorrected for bias errors and index E measured vari- 
able. Finally, x , y ,z and x ol y Q , z. are the position coordinates 
a a a 8 8 8 

of oc and 8 vane with respect to airplane center of gravity. 
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PROBLEM STATEMENT 


The general form of state equations for a given system can be written as 
x(t) = £ [x(t), g(t) , © 1 ] + g [x( t) ] £(t) , x(0) = x 0 (5) 


and the discrete form of the measured equations as 

z(i) = h [x(i),Ti(i), 0^ + n ( i ) , i = 1 , 2, ...» N (6) 

where x, g, and z are the state, input and measurement vector 

respectively, is the vector of unknown biases and scale factor errors, Xq 

is the vector of unknown initial conditions, E, and n are the process and 
measurement noise vectors respectively, and N is the number of data points. 

The compatibility check can be now formulated as an identification problem 
which involves the estimation of state and output variables, unknown para- 
meters 0 1 and Xq, and covariance matrices of £ and n, from measured data. 

The postulated model equations represent a nonlinear stochastic system 
with state-dependent process noise and with nonlinear output equations with an 
additive measurement noise. The state estimation in this case would be an 
extremely difficult problem. The separate estimation of unknown parameters 
would be equally complicated because of the resulting form of the sensitivity 
equations. For these reasons, possible simplification of the problem will be 
considered . 


ESTIMATION METHOD 

The state and parameter estimation problem outlined above can be reduced 
to parameter estimation only by neglecting the process noise altogether. The 
state estimation is thus replaced by the integration of the state equations 
with the estimated values for initial conditions and bias errors. For the 
parameter estimation the maximum likelihood (ML) method is applied. The 
measurement noise in equation (6) is assumed to be zero-mean, uncorrelated, 
and gaussian, i.e., 

E(n(i)} = 0, and E{n(i)n T (j)} = R(g 

where the symbol 6. . is the Kronecker delta. 

i/ 3 

The ML method finds a set of parameters by minimizing the log - 
likelihood function (see, e.g. reference 14). 

v T (i) R -1 v(i) - | Jin | R | (7) 
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where 


» - h *J T 

v(i ) = z(i) - h [x(i ) f ri(i ) » 0 ] 


( 8 ) 


2 

a is the variance of the measurement noise in output variables and N is 

y 

number of data points. Minimizing (7) for parameters in R gives the estimate 
of measurement-noise covariance matrix as 


1 


N 


E = “ if, 


v(i ) v (i ) 


(9) 


The estimates of the remaining unknown parameters are given by the root of the 
equation 

3J ( 0) I 


3 0 


= 0 


( 10 ) 


0=0 


for R replaced by R. This root can be found by modified Newton-Raphson itera- 
tion technique (ref. 14). 


It is well known that under the above mentioned assumptions the final 
estimates of unknown parameters are consistent, asymptotically unbiased, and 
asymptotically efficient with 


E { ( 0 - 0) (0 - 0) T } > Z [a T R~V] (11) 

i=1 

T 

In this expression A is the transpose of the sensitivity matrix with the 
elements 





1 , 2 , 


m i 


1 , 2 , . . . . n 


where n m and np are the numbers of measured output variables and unknown 
parameters in 0 respectively. 


COMPUTING ALGORITHM 

The block diagram of computing procedure for the ML estimation of unknown 
parameters is presented in Figure 1 . The measured data are considered in the 
form of digitized time histories of input and output variables with sampling 
interval At. The state and output equations are given by Egs. (2) and (3) 
respectively with the unknown biases and scale factors 
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of= [b ,b , b ,b , b , b ,b ,b ,b ,b.,b ,A , 
1 L ax ay az p q r v g a <J> 9 v 

VVV VV V 


and unknown initial conditions 


x o = tvvvv v V’U 

For the integration of state equations and computing output time histories the 
starting values of unknown parameters must be specified. The biases and scale 
factors are usually set equal to zero. The starting value of initial condi- 
tions are obtained from measured variables at t=0 as 


U 0 = V R,E (0) COS a R,E ( ° ) COSf3 R,E (0) 


V 0= V R,E (0) Sin0 R,E (O) 

V R,E (0) Sina R f E (0) COS0 R >E (O) 

( 12 ) 

h R,E (0) 

V *R,E (0) 


V Ve <0) 

Using the modified Newton -Raphson iterative technique the unknown para- 
meters obtained from 

0 = G X1 + AO (13) 

N 



where 0., are the starting values 

and N 

N -1 N 

40 = [ E A. R A. 1 f E A. R V. 

.,1 1 JL .,1 

1=1 1=1 

The transpose of the sensitivity matrix A is given in Table I, where the 
crosses indicate sensitivities 3y , /30 = -3v /30 computed by a numerical 

X* jC X 

method of reference 14. The remaining sensitivities in Table I are equal to a 
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known constant or zero, or are formed by the computed output variables 
estimates of elements in the measurement noise covariance matrix R are 
obtained from Eg. (9) and the residuals from Eg. (8). 


The 


For t=0 the output eguations for <)> and 9 give 

R R 

V *R (0) - (1+ V *0 

b 9 = e R ( o ) - (i + x 9 )e 0 


(15) 


It means that for X. and X a a linear relationship exists between 

<() a 

and 0 and b a . The algorithm therefore provides two options for 
0 0 

being either estimated or computed from Eg. (15). 


♦o and V 

b, and b Q 
(() 0 


The iteration process is completed when certain stopping criterion is 
met. In this report two criteria proposed in reference 13 are adopted, i.e.. 


A0 
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AJ 
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where and are specified, for example e 2 = .001. These two criteria 

should be met simultaneously. When the iteration is completed the correlation 
matrix of unknown parameters is computed from 


where 


M*L = WM 1 W 
MOD 


M = {m } 
w = {1/^77} 


(16) 


The residuals can be examined using their time histories, autocorrelation 
functions and power spectral densities (ref. 9). 

EXAMPLES USING SIMULATED DATA 

The method developed in this study was applied to simulated data to check 
the accuracy of the computing algorithm and estimated parameters. Two sets of 
simulated data, one representing the longitudinal and the other representing 
the lateral motion of an airplane, were used. The time histories of the input 
and output variables are shown in Figures 2 and 3. The sampling interval for 
these time histories was At= 0.05 sec. The measurement noise with standard 
errors given in Table II was added to the simulated data. 
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The results of three cases using the longitudinal data with measurement 
noise only in the output variables are presented in Table III. In Case 1 the 
initial conditions were fixed on true values, in Case 2 on values corrupted by 
measurement noise. In Case 3 the initial conditions were estimated. The 
degradation in the accuracy of the estimated parameters from Case 1 to Cases 2 
and 3 is apparent. It is probably caused by incorrect values of initial 
conditions (Case 2), by an increased number of unknown parameters (Case 3) and 
high correlation (greater than 0.9) between some parameters. The expected 
high correlation between parameters b^ and 0^ (see Eg. 15) did not materialize. 

A similar approach was adopted in the next three cases, where both the 
input and output variables were corrupted by measurement noise and the number 
of unknown parameters was increased. The results in Table IV do not show any 
substantial deterioration in the accuracy of the estimated parameters when 
compared with the previous results in Case 2 and 3. 

In Table V two sets of results from the analysis of lateral data are 
given. These sets differ in the two values of simulated bias errors in the 
variables 0, p, r, and <j> . The different values of bias errors did not 
change the differences between the true and mean values of estimated para- 
meters but changed the standard errors of the parameters and the correlation 
between them. The accuracy was worse than that obtained from the longitudinal 
data. The estimation was then repeated with the initial conditions as the 
unknown parameters. This attempt, however, failed because of the estimation 
procedure divergence. The analysis of lateral data indicated a possibility of 
identif iability problems with this type of data. To avoid that, the design of 
an optimal maneuver for more accurate estimates and with less sensitivity to 
the number of unknown parameters should be investigated. 


ANALYSIS OF FLIGHT DATA 

Three sets of measured flight data were analyzed. Two of them repre- 
sented a longitudinal motion of an airplane, the last set was obtained from a 
combined maneuver with predominantly lateral motion. The sampling interval 
for all data was At= 0.05 sec. The measured output variables V R , , 

and 0 in the first run are presented in Figure 4. The resulting estimates 
R 

which include the parameter mean values, their standard errors (Cramer-Rao 
lower bound) and standard errors of the measurement noise in the output vari- 
ables are summarized in Table VI. In Case 1 the vector of unknown parameters 
was postulated as 



,b 


az 


b 

Q 



V 


X ,u 
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0 


,w 


0 ' 



As can be seen from the results, three pairs of estimated parameters are 
highly correlated, and the parameters b and 0 , have large standard errors. 

d a (J 

As the next step, therefore, the parameters b ax , Wq and X were fixed on their 
estimated values. The estimation of the remaining parameters was then repeated 
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in Case 2. The new results indicate no change in the mean values but lower 
standard errors of the estimates. The predicted time histories of the output 
variables are compared with those measured in Figure 4. The agreement is very 
good in all variables plotted. 

The time histories in the second run had similar form as those in the 
first one. The estimated parameters are contained in the last two columns of 
Table VI. All parameters from both runs agree well with the exception of the 
parameter by. This disagreement could be caused by the elimination of Ay in 

the vector of unknown parameters in the second run. As in Run 1 the results 
also exhibit high correlation between some parameters. 

In order to further validate the results, the ML estimates from Run 1 
were compared in Table VII with those obtained by a nonlinear-fixed-internal 
smoothing (NFIS) technique (ref. 12) and by an extended Kalman filter (EKF) 
(ref. 9). For the ML and NFIS estimation the initial values of unknown bias 
and scale factor errors were set equal to zero whereas for the EKF these 
values were made equal to the ML estimates because of the slow convergence of 
the filter. The estimates from all the three techniques aqree well. The main 
differences are seen only in parameters by and 0^. The parameter 9 q was 

estimated with poor accuracy in all cases. The reason for the disagreement in 

by could be due to insufficient excitation of the airspeed during the airplane 

motion. 

Data from the third maneuver were analyzed assuming only the bias errors 
in variables a y , a z , p, q, r, 0, <)> and 0. After reviewing the estimates, some 

of the less important terms were eliminated. As a result of that, the vector 
of unknown parameters was postulated as 



b 

az 


WV u o 


vv 



with b and b Q computed from Eq. (15). The measured and predicted time 
cp a 

histories of output variables are plotted in Figure 5. Then, as in the pre- 
vious case, the ML estimates were compared with those using NFIS and EKF 
techniques. The results are given in Table VIII. The agreement between the 
ML and EKF estimates is good. This could be, however, due to the use of ML 
estimates as starting values for the EKF technique. Therefore, more thorough 
checks should be made for better assessment of results from both methods. The 
main differences between ML and NFIS techniques are in the parameters b^, and 

initial conditions v 0 ,<|> 0 and 0 Q . These differences can be caused by poor 

accuracy of these estimates and different values of the remaining parameters. 
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CONCLUDING REMARKS 


A maximum likelihood method was developed for the estimation of initial 
conditions and bias errors in measured airplane responses. In the development 
it was assumed that the input variables to the system represented by linear 
accelerations and angular velocities are measured without random errors. The 
model relating airplane state, input and output variables is based on six- 
degree -of -freedom kinematic equations and on output equations specifying the 
measured variables. 

The resulting technique was first applied to a limited number of simu- 
lated data runs to check the accuracy of the computing algorithm and estimated 
parameters. It was demonstrated that the increased number of unknown para- 
meters and the correlation among them can degrade the accuracy of the esti- 
mates. At the same time it was observed that a moderate noise in measured 
inputs has only a small effect on the accuracy of the results. The lateral 
maneuver analyzed provided less accurate parameter estimates than the longi- 
tudinal one. The increased number of unknowns in the lateral case resulted in 
a divergence of the estimation procedure. The results from simulated data 
indicate a need for a design of an optimal maneuver for more accurate esti- 
mates and with less sensitivity to the number of unknown parameters. 

The maximum likelihood method developed for this report was also applied 
to the analysis of real flight data. Two longitudinal maneuvers similar in 
form were analyzed first. The resulting parameter estimates from both maneu- 
vers were in good agreement. Then the maximum likelihood estimates were 
compared with those obtained by a nonlinear-fixed-interval smoother and an 
extended Kalman filter. The comparison of the three techniques showed no main 
differences in the results. Similar conclusions were obtained from the anal- 
ysis of a lateral maneuver with a strong longitudinal coupling. All the 
comparisons served as a verification of the maximum likelihood technique 
presented in the report. 
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TABLE I. - TRANSPOSE OF SENSITIVITY MATRIX 



x indicate sensitivities computed by a numerical method 
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TABLE II. - STANDARD ERRORS OF 
MEASUREMENT NOISE 


Standard Error 
of Measurement 

Variable Noise of Variable 


a x , m/sec^ 0.1 

o 

a y/ m/sec 0.1 

a 2 , m/sec 0.1 

p, rad/sec 0.002 

q, rad/sec 0.002 

r, rad/sec 0.002 

V, m/s ec 1 • 0 

3, rad 0.002 

а, rad 0.002 

rad 0.002 

б, rad 0.002 

\p, rad 0.002 
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TABLE III. - EFFECT OF MEASUREMENT NOISE IN OUTPUT VARIABLES ON 
PARAMETER ESTIMATES FOR SIMULATED DATA OF LONGITUDINAL MOTION 
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with high correlation 



TABLE IV. - EFFECT OF MEASUREMENT NOISE IN INPUT AND OUTPUT VARIABLES 
ON PARAMETER ESTIMATES FOR SIMULATED DATA OF LONGITUDINAL MOTION 



nj 

Vj Ci 
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with high correlation 



TABLE V. - EFFECT OF MEASUREMENT NOISE ON PARAMETER 
ESTIMATES FOR SIMULATED DATA OF LATERAL MOTION 



17 








TABLE VI. - ESTIMATES OF PARAMETERS FROM LONGITUDINAL 
FLIGHT DATA OF TWO DIFFERENT RUNS 



a Parameters with high correlation 
k Computed from 9^ and 9p^g(0) 
c Fixed parameter 


I 
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TABLE VII. - ESTIMATES OF PARAMETERS FROM LONGITUDINAL 
FLIGHT DATA USING THREE DIFFERENT METHODS 


Parameter 

ML 

NFIS 

EKF 

0 

s { 0) 

A 

0 

s ( 0) 

0 

s (0) 

b a2 / 9 units 

-.1240 


.00021 

-.122 

.0035 

-.123 

.0035 

bg, deg/sec 

1 .792 

a 

.0018 

1 .8 

.14 

1 .74 

.042 

by, m/sec 

1 .24 

a 

.098 

.02 

.40 

-.28 

.075 

b , deg 
a 

2.631 


.0036 

2.34 

.34 

2.6 

.10 

bg, deg 

-12.72 

b 


-13.03 

.10 

-12.64 

.033 

\ 

-.073 


.0026 

-.050 

.015 

-.09 

.012 

a 








Uq, m/sec 

35.5 

a 

.10 

37.8 

.40 

35.5 

.82 

w 0 , m/sec 

6.0 

c 


6.1 

.27 

6.1 

.39 

deg 

-.19 

a 

.025 

.7 

.30 

.4 

.72 

s(V), m/sec 

.260 



.177 




s ( a) , deg 

.280 



.158 




s ( 6 ) , deg 

.748 



1 .08 





a Parameters with high correlation 
13 Computed from 0 and 8 (0) 

r» . 0 K / r« 

Fixed parameter 
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TABLE VIII. - ESTIMATES OF PARAMETERS FROM 
FLIGHT DATA USING THREE DIFFERENT METHODS 



a Parameter with high correlation 
b Computed from <j> Qf 8 Q and ♦ R ,E (0) ' 9 r,e (0) 
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n R,E 


I L_ 

Figure 1. 



- Block diagram of computing procedure for 
ML estimation of unknown parameters. 
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Figure 3. - Simulated lateral motion of an airplane 
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